library(readxl)
library("minpack.lm", lib.loc="/Library/Frameworks/R.framework/Versions/3.4/Resources/library")
randomisation_wavelength_Autosaved_ <- read_excel("/Users/quentinguignard/Desktop/Chapter 3 Sirex eyes/wurzburg/randomisation wavelength (Autosaved).xlsx", 
                                                  sheet = "data_ERG_R")

#####################################################
###### Data mangaement 
#####################################################
Dark = randomisation_wavelength_Autosaved_[1:17,]
Dim = randomisation_wavelength_Autosaved_[18:34,-c(3,4,13,14,15)]
Strong = randomisation_wavelength_Autosaved_[35:51, -c(3,4,11,13,14,15)]


Dark_females = Dark[,1:12]
Dark_males = Dark[,c(1,2,13:22)]
Dim_females = Dim[,1:10]
Dim_males = Dim[,c(1,2,11:17)]
Strong_females = Strong[,1:9]
Stong_males = Strong[,c(1,2,10:16)]
Dark_alpha_green = Dark[7:17,]
Strong_alpha_UV = Strong[1:7,]

mean_both_dark_alpha_green = data.frame(Dark[,2],apply(Dark[7:17,3:22], 1, mean), apply(Dark[7:17,3:22], 1, sd))
colnames(mean_both_dark_alpha_green) = c('Wavelength', "mean_both", "sd_both")
mean_male_dark_alpha_green = data.frame(Dark[7:17,2],apply(Dark[7:17,c(13:22)], 1, mean), apply(Dark[7:17,c(13:22)], 1, sd))
colnames(mean_male_dark_alpha_green) = c('Wavelength', "mean_males", "sd_males")
mean_female_dark_alpha_green = data.frame(Dark[7:17,2],apply(Dark[7:17,c(3:12)], 1, mean), apply(Dark[7:17,c(3:12)], 1, sd))
colnames(mean_female_dark_alpha_green) = c('Wavelength', "mean_females", "sd_females")

mean_both_Strong_alpha_UV = data.frame(Strong_alpha_UV[,2],apply(Strong_alpha_UV[,3:16], 1, mean), apply(Strong_alpha_UV[,3:16], 1, sd))
colnames(mean_both_Strong_alpha_UV) = c('Wavelength', "mean_both", "sd_both")
mean_male_Strong_alpha_UV = data.frame(Strong_alpha_UV[,2],apply(Strong_alpha_UV[,10:16], 1, mean), apply(Strong_alpha_UV[,10:16], 1, sd))
colnames(mean_male_Strong_alpha_UV) = c('Wavelength', "mean_males", "sd_males")
mean_female_Strong_alpha_UV = data.frame(Strong_alpha_UV[,2],apply(Strong_alpha_UV[,3:9], 1, mean), apply(Strong_alpha_UV[,3:9], 1, sd))
colnames(mean_female_Strong_alpha_UV) = c('Wavelength', "mean_females", "sd_females")

mean_both_dark = data.frame(Dark[,2],apply(Dark[,3:22], 1, mean), apply(Dark[,3:22], 1, sd))
colnames(mean_both_dark) = c('Wavelength', "mean_both", "sd_both")
mean_both_Strong = data.frame(Strong[,2],apply(Strong[,3:16], 1, mean), apply(Strong[,3:16], 1, sd))
colnames(mean_both_Strong) = c('Wavelength', "mean_both", "sd_both")
mean_both_Dim = data.frame(Dim[,2],apply(Dim[,3:17], 1, mean), apply(Dim[,3:17], 1, sd))
colnames(mean_both_Dim) = c('Wavelength', "mean_both", "sd_both")

mean_male_Strong_UV = data.frame(Stong_males[,2],apply(Stong_males[,3:9], 1, mean), apply(Stong_males[,3:9], 1, sd))
colnames(mean_male_Strong_UV) = c('Wavelength', "mean_males", "sd_males")
mean_female_Strong_UV = data.frame(Strong_females[,2],apply(Strong_females[,3:9], 1, mean), apply(Strong_females[,3:9], 1, sd))
colnames(mean_female_Strong_UV) = c('Wavelength', "mean_females", "sd_females")

Dark_mean_males = data.frame(Dark[,2],apply(Dark_males[,3:12], 1, mean), apply(Dark_males[,3:12], 1, sd))
colnames(Dark_mean_males) = c('Wavelength', "mean_males", "sd_both")
Dark_mean_females = data.frame(Dark[,2],apply(Dark_females[,3:12], 1, mean), apply(Dark_females[,3:12], 1, sd))
colnames(Dark_mean_females) = c('Wavelength', "mean_females", "sd_both")

Dim_mean_males = data.frame(Dim[,2],apply(Dim_males[,3:9], 1, mean), apply(Dim_males[,3:9], 1, sd))
colnames(Dim_mean_males) = c('Wavelength', "mean_males", "sd_both")
Dim_mean_females = data.frame(Dim[,2],apply(Dim_females[,3:9], 1, mean), apply(Dim_females[,3:9], 1, sd))
colnames(Dim_mean_females) = c('Wavelength', "mean_males", "sd_both")

#####################################################
###### Alpha peak of green peak 
#####################################################

model_green_alpha_male_all_spectrum_2parameter_fix = nlsLM(mean_males ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                           start = list( x=575), data = Dark_mean_males)
model_green_alpha_female_all_spectrum_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                             start = list(x=575), data = Dark_mean_females)

model_green_alpha_male_right_side2_2parameter_fix = nlsLM(mean_males ~ exp((-396*(log10(Wavelength/x))^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                          start = list(x=575), data = Dark_mean_males[13:17,])
model_green_alpha_female_right_side2_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                            start = list( x=575), data = Dark_mean_females[13:17,])

model_green_alpha_male_alpha_green_2parameter_fix = nlsLM(mean_males ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                          start = list( x=575), data = Dark_mean_males[7:17,])
model_green_alpha_female_alpha_green_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                            start = list( x=575), data = Dark_mean_females[7:17,])



model_green_alpha_male_all_spectrum_4parameter_fix = nlsLM(mean_males ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                           start = list( x=575), data = Dark_mean_males)
model_green_alpha_female_all_spectrum_4parameter_fix = nlsLM(mean_females ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                             start = list(x=575), data = Dark_mean_females)

model_green_alpha_male_right_side2_4parameter_fix = nlsLM(mean_males ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                          start = list(x=575), data = Dark_mean_males[13:17,])
model_green_alpha_female_right_side2_4parameter_fix = nlsLM(mean_females ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                            start = list( x=575), data = Dark_mean_females[13:17,])

model_green_alpha_male_alpha_green_4parameter_fix = nlsLM(mean_males ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                          start = list( x=575), data = Dark_mean_males[7:17,])
model_green_alpha_female_alpha_green_4parameter_fix = nlsLM(mean_females ~ exp((-387*log10(Wavelength/x)^2)*(1+14.0*log10(Wavelength/x) + 120*(log10(Wavelength/x))^2 + 418*(log10(Wavelength/x))^3)), 
                                                            start = list( x=575), data = Dark_mean_females[7:17,])


model_green_alpha_male_all_spectrum_Govardovskii = nlsLM(mean_males ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674),
                                                           start = list( x=550), data = Dark_mean_males)
model_green_alpha_female_all_spectrum_Govardovskii = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                             start = list(x=575), data = Dark_mean_females)

model_green_alpha_male_right_side2_Govardovskii = nlsLM(mean_males ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                          start = list(x=575), data = Dark_mean_males[13:17,])
model_green_alpha_female_right_side2_Govardovskii = nlsLM(mean_females ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                            start = list( x=575), data = Dark_mean_females[13:17,])

model_green_alpha_male_alpha_green_Govardovskii = nlsLM(mean_males ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                          start = list( x=575), data = Dark_mean_males[7:17,])
model_green_alpha_female_alpha_green_Govardovskii = nlsLM(mean_females ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                            start = list( x=575), data = Dark_mean_females[7:17,])

summary(model_green_alpha_male_all_spectrum_2parameter_fix)
summary(model_green_alpha_female_all_spectrum_2parameter_fix)
summary(model_green_alpha_male_right_side2_2parameter_fix)
summary(model_green_alpha_female_right_side2_2parameter_fix)
summary(model_green_alpha_male_alpha_green_2parameter_fix)
summary(model_green_alpha_female_alpha_green_2parameter_fix)

summary(model_green_alpha_male_all_spectrum_Govardovskii)
summary(model_green_alpha_female_all_spectrum_Govardovskii)
summary(model_green_alpha_male_right_side2_Govardovskii)
summary(model_green_alpha_female_right_side2_Govardovskii)
summary(model_green_alpha_male_alpha_green_Govardovskii)
summary(model_green_alpha_female_alpha_green_Govardovskii)


predict.Wavelengh[which.max(predict(model_green_alpha_male_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_green_alpha_male_alpha_green_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_alpha_green_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]


predict.Wavelengh[which.max(predict(model_green_alpha_male_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_green_alpha_male_right_side2_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_right_side2_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_green_alpha_male_alpha_green_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_green_alpha_female_alpha_green_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]

plot(Dark_mean_males[,2] ~ Dark_mean_males$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green1")
lines(330:660, predict(model_green_alpha_male_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green2")
lines(330:660, predict(model_green_alpha_male_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(330:660))), col = "green3")
lines(330:660, predict(model_green_alpha_male_all_spectrum_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green3")


plot(Dark_mean_females[,2] ~ Dark_mean_females$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_green_alpha_female_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green1")
lines(330:660, predict(model_green_alpha_female_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green2")
lines(330:660, predict(model_green_alpha_female_alpha_green_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "green3")


plot(Dark_mean_males[,2] ~ Dark_mean_males$Wavelength, type = 'l', ylim=c(0,1))
lines(Dark_mean_males$Wavelength, predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(Dark_mean_males$Wavelength))))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_green_alpha_male_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_green_alpha_male_alpha_green_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))      
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")


plot(Dark_mean_females[,2] ~ Dark_mean_females$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_green_alpha_female_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_green_alpha_female_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_green_alpha_female_alpha_green_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))      
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")


#####################################################
###### UV peak
#####################################################

model_UV_alpha_male_all_spectrum_2parameter_fix = nlsLM(mean_males ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                           start = list( x=350), data = mean_male_Strong_UV)

model_UV_alpha_female_all_spectrum_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                             start = list(x=350), data = mean_female_Strong_UV)

model_UV_alpha_male_alpha_UV_2parameter_fix = nlsLM(mean_males ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                          start = list( x=350), data = mean_male_Strong_UV[1:7,])

model_UV_alpha_female_alpha_UV_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                            start = list( x=350), data = mean_female_Strong_UV[1:7,])

model_UV_alpha_male_left_side_UV_2parameter_fix = nlsLM(mean_males ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                    start = list( x=350), data = mean_male_Strong_UV[1:4,])

model_UV_alpha_female_left_side_UV_2parameter_fix = nlsLM(mean_females ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                      start = list( x=350), data = mean_female_Strong_UV[1:4,])



model_UV_alpha_male_all_spectrum_Govardovskii = nlsLM(mean_males ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                        start = list( x=350), data = mean_male_Strong_UV)

model_UV_alpha_female_all_spectrum_Govardovskii = nlsLM(mean_females ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                          start = list( x=350), data = mean_female_Strong_UV)

model_UV_alpha_male_alpha_UV_Govardovskii = nlsLM(mean_males ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                    start = list( x=350), data = mean_male_Strong_UV[1:7,])

model_UV_alpha_female_alpha_UV_Govardovskii = nlsLM(mean_females ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674), 
                                                      start = list( x=350), data = mean_female_Strong_UV[1:7,])

plot(mean_male_Strong_UV[,2] ~ mean_male_Strong_UV$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_UV_alpha_male_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple1")
lines(330:660, predict(model_UV_alpha_male_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple3")
lines(330:660, predict(model_UV_alpha_male_all_spectrum_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple", lty = 2)
predict.Wavelengh = c(330:660)

predict.Wavelengh[which.max(predict(model_UV_alpha_male_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_UV_alpha_male_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_UV_alpha_male_left_side_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_left_side_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]


predict.Wavelengh[which.max(predict(model_UV_alpha_male_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_UV_alpha_male_alpha_UV_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_alpha_UV_Govardovskii, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_UV_alpha_male_all_spectrum_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_all_spectrum_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

predict.Wavelengh[which.max(predict(model_UV_alpha_male_alpha_UV_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]
predict.Wavelengh[which.max(predict(model_UV_alpha_female_alpha_UV_4parameter_fix, newdata = data.frame(Wavelength = c(330:660))))]

summary(model_UV_alpha_male_all_spectrum_2parameter_fix)
summary(model_UV_alpha_female_all_spectrum_2parameter_fix)
summary(model_UV_alpha_male_alpha_UV_2parameter_fix)
summary(model_UV_alpha_female_alpha_UV_2parameter_fix)
summary(model_UV_alpha_male_left_side_UV_2parameter_fix)
summary(model_UV_alpha_female_left_side_UV_2parameter_fix)

summary(model_UV_alpha_male_all_spectrum_Govardovskii)
summary(model_UV_alpha_female_all_spectrum_Govardovskii)
summary(model_UV_alpha_male_alpha_UV_Govardovskii)
summary(model_UV_alpha_female_alpha_UV_Govardovskii)

summary(model_UV_alpha_male_all_spectrum_4parameter_fix)
summary(model_UV_alpha_female_all_spectrum_4parameter_fix)
summary(model_UV_alpha_male_alpha_UV_4parameter_fix)
summary(model_UV_alpha_female_alpha_UV_4parameter_fix)

plot(mean_female_Strong_UV[,2] ~ mean_female_Strong_UV$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_UV_alpha_female_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple")
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple")
lines(330:660, predict(model_UV_alpha_male_left_side_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple")

#####################################################
###### graph opsins
#####################################################

plot(Dark_mean_males[,2] ~ Dark_mean_males$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_UV_alpha_male_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple")


plot(Dark_mean_females[,2] ~ Dark_mean_females$Wavelength, type = 'l', ylim=c(0,1))
lines(330:660, predict(model_green_alpha_female_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "green1")
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple")

quartz(width = 13, height = 8, bg= 'white')
plot(mean_male_Strong_UV[,2] ~ mean_male_Strong_UV$Wavelength, lwd = 2, cex.lab=1.5, type = 'l', ylim=c(0,1.1), col = 'blue', xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"))
arrows(mean_male_Strong_UV$Wavelength, mean_male_Strong_UV[,2]-mean_male_Strong_UV[,3], mean_male_Strong_UV$Wavelength, mean_male_Strong_UV[,2]+mean_male_Strong_UV[,3],
       length=0.05, angle=90, code=3, col = 'blue', lwd = 2)
lines(mean_female_Strong_UV[,2] ~ mean_female_Strong_UV$Wavelength, lwd = 2, type = 'l', ylim=c(0,1), col = 'pink')
arrows(mean_female_Strong_UV$Wavelength, mean_female_Strong_UV[,2]-mean_female_Strong_UV[,3], mean_female_Strong_UV$Wavelength, mean_female_Strong_UV[,2]+mean_female_Strong_UV[,3],
       length=0.05, angle=90, code=3, col = 'pink', lwd = 2)
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "purple", lty = 2)
title(main = "Relative sensitivity (±Sd) of Males (n=7) and Females (n=7) after strong green light adaptation",cex.main=1.5)
legend("topright", legend = c('Males', 'Females', 'UV opsin (max = 363)'),cex=1.5, lty = c(1,1,2), lwd = 2, col = c('blue', 'pink', "purple"))

quartz(width = 13, height = 8, bg= 'white')
plot(Dark_mean_males[,2] ~ Dark_mean_males$Wavelength, type = 'l', cex.lab=1.5, ylim=c(0,1.1), col = 'blue', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), lwd = 2)
arrows(Dark_mean_males$Wavelength, Dark_mean_males[,2]-Dark_mean_males[,3], Dark_mean_males$Wavelength, Dark_mean_males[,2]+Dark_mean_males[,3],
       length=0.05, angle=90, code=3, col = 'blue', lwd = 2)
lines(Dark_mean_females[,2] ~ Dark_mean_females$Wavelength, type = 'l', ylim=c(0,1), col = 'pink2', lwd = 2)
arrows(Dark_mean_females$Wavelength, Dark_mean_females[,2]-Dark_mean_females[,3], Dark_mean_females$Wavelength, Dark_mean_females[,2]+Dark_mean_females[,3],
       length=0.05, angle=90, code=3, col = 'pink2', lwd = 2)
lines(330:660, predict(model_green_alpha_male_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "Chartreuse 1", lty = 2, lwd = 2)
lines(330:660, predict(model_green_alpha_female_right_side2_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "Chartreuse 4", lty = 4, lwd = 2)
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple", lwd = 2, lty = 2)
title(main = "Relative sensitivity (±Sd) of Males (n=10) and Females (n=10) after dark adaptation",cex.main=1.5)
legend("topright", legend = c('Males', 'Females', 'UV opsin (max = 363)', 'Green opsin males (max = 520)','Green opsin females (max = 527)' ),
       lty = c(1,1,2,2,4), col = c('blue', 'pink','purple', "Chartreuse 1", 'Chartreuse 4'), lwd = 2,cex=1.10)


quartz(width = 13, height = 8, bg= 'white')
plot(mean_male_Strong_UV[,2] ~ mean_male_Strong_UV$Wavelength, lwd = 2, cex.lab=1.5, type = 'l', ylim=c(0,1.1), col = 'blue', xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"))
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "purple1", lty = 'dashed')
lines(320:650, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "purple2", lty = 'dotted')
lines(340:670, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "purple3", lty = 'twodash')







####################

quartz(width = 13, height = 8, bg= 'white')
plot(Dim_mean_males[,2] ~ Dim_mean_males$Wavelength, type = 'l', cex.lab=1.5, ylim=c(0,1.1), col = 'blue', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), lwd = 2)
arrows(Dim_mean_males$Wavelength, Dim_mean_males[,2]-Dim_mean_males[,3], Dim_mean_males$Wavelength, Dim_mean_males[,2]+Dim_mean_males[,3],
       length=0.05, angle=90, code=3, col = 'blue', lwd = 2)
lines(Dim_mean_females[,2] ~ Dim_mean_females$Wavelength, type = 'l', ylim=c(0,1), col = 'pink2', lwd = 2)
arrows(Dim_mean_females$Wavelength, Dim_mean_females[,2]-Dim_mean_females[,3], Dim_mean_females$Wavelength, Dim_mean_females[,2]+Dim_mean_females[,3],
       length=0.05, angle=90, code=3, col = 'pink2', lwd = 2)
lines(330:660, predict(model_green_alpha_male_right_side2_2parameter_fix_dim_light, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "Chartreuse 1", lty = 2, lwd = 2)
lines(330:660, predict(model_green_alpha_female_right_side2_2parameter_fix_dim_light, newdata = data.frame(Wavelength = c(330:660)))
      + 0.29*exp((-247*log10(c(330:660)/350 )^2)*(1+3.59*log10(c(330:660)/350 ) + (3/8)*(3.59*log10(c(330:660)/350))^2))
      , col = "Chartreuse 4", lty = 4, lwd = 2)
lines(330:660, predict(model_UV_alpha_female_alpha_UV_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), col = "purple", lwd = 2, lty = 2)
title(main = "Relative sensitivity (±Sd) of Males (n=7) and Females (n=7) after dim green adaptation",cex.main=1.5)
legend("topright", legend = c('Males', 'Females', 'UV opsin (max = 363)', 'Green opsin males (max = 520)','Green opsin females (max = 527)' ),
       lty = c(1,1,2,2,4), col = c('blue', 'pink','purple', "Chartreuse 1", 'Chartreuse 4'), lwd = 2,cex=1.10)

model_green_alpha_male_right_side2_2parameter_fix_dim_light = nlsLM(mean_males ~ exp((-380*(log10(Wavelength/x))^2)*(1+6.09*log10(Wavelength/x) + (3/8)*(6.09*log10(Wavelength/x))^2)), 
                                                          start = list(x=575), data = Dark_mean_males[13:17,])

model_green_alpha_female_right_side2_2parameter_fix_dim_light = nlsLM(mean_females ~ exp((-380*log10(Wavelength/x)^2)*(1+6.09*log10(Wavelength/x) + (3/8)*(6.09*log10(Wavelength/x))^2)), 
                                                            start = list( x=575), data = Dark_mean_females[13:17,])
abline(model_green_alpha_male_right_side2_2parameter_fix_dim_light)


#####################################################
###### graph opsins fitting
#####################################################
Wavelength = c(334, 353, 361, 382, 399, 421, 432, 453, 471, 504, 521, 532, 554, 572, 603, 623, 652)
Intensity = c(0.936, 1.03, 0.95, 1.12, 1.11, 1.07, 1.06, 1.11, 1.15, 1.09, 0.974, 1.05, 1.14, 0.973, 1.05, 1.05, 1.01)

model_green_alpha_both_all_spectrum_2parameter_fix = nlsLM(mean_both ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                           start = list( x=575), data = mean_both_dark[13:17,])
model_green_alpha_both_all_spectrum_Govardovskii = nlsLM(mean_both ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674),
                                                         start = list( x=550), data = mean_both_dark[13:17,])
model_UV_alpha_both_all_spectrum_2parameter_fix = nlsLM(mean_both ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                        start = list( x=350), data = mean_both_Strong)
model_UV_alpha_both_all_spectrum_Govardovskii = nlsLM(mean_both ~ 1/(exp(69.7*(0.8892512-(Wavelength/x)^-1)) + exp(28*(0.922-(Wavelength/x)^-1)) + exp(-14.9*(1.104-(Wavelength/x)^-1)) + 0.674),
                                                         start = list( x=350), data = mean_both_Strong)
Value_waveleength_tested_model_green_alpha_both_all_spectrum_2parameter_fix = c(predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(334, 353, 361, 382, 399, 421, 432, 453, 471, 504, 521, 532, 554, 572, 603, 623, 652))))
Value_waveleength_tested_model_Govardovskii_green_alpha_both_all_spectrum_2parameter_fix = c(predict(model_green_alpha_both_all_spectrum_Govardovskii, newdata = data.frame(Wavelength = c(334, 353, 361, 382, 399, 421, 432, 453, 471, 504, 521, 532, 554, 572, 603, 623, 652))))
Value_waveleength_tested_model_UV_alpha_both_all_spectrum_2parameter_fix = c(predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(334, 353, 361, 382, 399, 421, 432, 453, 471, 504, 521, 532, 554, 572, 603, 623, 652))))

model2_green_alpha_both_all_spectrum_2parameter_fix = nlsLM(mean_both ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                           start = list( x=575), data = Dark[14:17,3:22])
model2_UV_alpha_both_all_spectrum_2parameter_fix = nlsLM(mean_both ~ exp((-396*log10(Wavelength/x)^2)*(1+10.1*log10(Wavelength/x) + (3/8)*(10.1*log10(Wavelength/x))^2)), 
                                                        start = list( x=350), data = mean_both_Strong)

scaling_factor_UV_dark = mean_both_dark$mean_both[which(mean_both_Strong$Wavelength == 361)]/ mean_both_Strong$mean_both[which(mean_both_dark$Wavelength == 361)]
scaling_factor_green_strong = mean_both_Strong$mean_both[which(mean_both_dark$Wavelength == 521)]/ mean_both_dark$mean_both[which(mean_both_dark$Wavelength == 521)]
scaling_factor_UV_dim = mean_both_Dim$mean_both[which(mean_both_Dim$Wavelength == 361)]/ mean_both_Strong$mean_both[which(mean_both_dark$Wavelength == 361)]
scaling_factor_green_dim = mean_both_Dim$mean_both[which(mean_both_Dim$Wavelength == 521)]/ mean_both_dark$mean_both[which(mean_both_dark$Wavelength == 521)]

quartz(width = 13, height = 8, bg= 'white')
par(mfrow = c(2,2), mar = c(5,4,4,2))
plot(mean_both_Strong[,2] ~ mean_both_Strong$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'grey', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_Strong$Wavelength, mean_both_Strong[,2]-mean_both_Strong[,3], mean_both_Strong$Wavelength, mean_both_Strong[,2]+mean_both_Strong[,3],
       length=0.05, angle=90, code=3, col = 'grey', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "purple", lty = "dashed")
lines(330:660,predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_green_strong, lwd = 2, col = "chartreuse3", lty = "dotted")
lines(330:660,predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_green_strong + predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) 
      , lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 14', cex = 1.1)
text(x = 327, y = 1.1, 'A', cex = 1.1, col = "black")
#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (7 males and 7 Females) 
#      after strong green light adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

#quartz(width = 13, height = 8, bg= 'white')
plot(mean_both_dark[,2] ~ mean_both_dark$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'grey', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_dark$Wavelength, mean_both_dark[,2]-mean_both_dark[,3], mean_both_dark$Wavelength, mean_both_dark[,2]+mean_both_dark[,3],
       length=0.05, angle=90, code=3, col = 'grey', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_UV_dark, lwd = 2, col = "purple", lty = 'dashed')
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "chartreuse3", lty = 'dotted')
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) + predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_UV_dark
      , lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 20', cex = 1.1)
text(x = 327, y = 1.1, 'B', cex = 1.1, col = "black")

#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (10 males and 10 Females) 
#      after dark adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

#quartz(width = 13, height = 8, bg= 'white')
plot(mean_both_Dim[,2] ~ mean_both_Dim$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'grey', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_Dim$Wavelength, (mean_both_Dim[,2]-mean_both_Dim[,3]), mean_both_Dim$Wavelength, (mean_both_Dim[,2]+mean_both_Dim[,3]),
       length=0.05, angle=90, code=3, col = 'grey', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_UV_dim , lwd = 2, col = "purple", lty = 'dashed')
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_green_dim , lwd = 2, col = "chartreuse3", lty = 'dotted')
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_UV_dim +  predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_green_dim, 
      lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 15', cex = 1.1)
text(x = 327, y = 1.1, 'C', cex = 1.1, col = "black")

#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (7 males and 8 Females) 
#      after dark adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

plot(x= c(330:660), predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "purple", lty = 'dashed', type = 'l',
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "chartreuse3", lty = 'dotted')
text(x = 610, y = 1, expression(paste(lambda, 'max = 363 ± 3')), cex = 1.1, col = 'purple')
text(x = 610, y = 0.92, expression(paste(lambda, 'max = 527 ± 2')), cex = 1.1, col = 'chartreuse3')
text(x = 327, y = 1.1, 'D', cex = 1.1, col = "black")

###### Black and white
#######################

quartz(width = 13, height = 8, bg= 'white')
par(mfrow = c(2,2), mar = c(5,4,4,2))
plot(mean_both_Strong[,2] ~ mean_both_Strong$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'black', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_Strong$Wavelength, mean_both_Strong[,2]-mean_both_Strong[,3], mean_both_Strong$Wavelength, mean_both_Strong[,2]+mean_both_Strong[,3],
       length=0.05, angle=90, code=3, col = 'black', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "black", lty = "dashed")
lines(330:660,predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_green_strong, lwd = 2, col = "black", lty = "dotted")
lines(330:660,predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_green_strong + predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) 
      , lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 14', cex = 1.1)
text(x = 327, y = 1.1, 'A', cex = 1.5, font = 2, col = "black")

#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (7 males and 7 Females) 
#      after strong green light adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

#quartz(width = 13, height = 8, bg= 'white')
plot(mean_both_dark[,2] ~ mean_both_dark$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'black', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_dark$Wavelength, mean_both_dark[,2]-mean_both_dark[,3], mean_both_dark$Wavelength, mean_both_dark[,2]+mean_both_dark[,3],
       length=0.05, angle=90, code=3, col = 'black', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_UV_dark, lwd = 2, col = "black", lty = "dashed")
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))), lwd = 2, col = "black", lty = "dotted")
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) + predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) * scaling_factor_UV_dark
      , lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 20', cex = 1.1)
text(x = 327, y = 1.1, 'B', cex = 1.5, font = 2)

#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (10 males and 10 Females) 
#      after dark adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

#quartz(width = 13, height = 8, bg= 'white')
plot(mean_both_Dim[,2] ~ mean_both_Dim$Wavelength, lwd = 2, cex.lab=1.5, type = 'p', ylim=c(0,1.1), col = 'black', 
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
arrows(mean_both_Dim$Wavelength, (mean_both_Dim[,2]-mean_both_Dim[,3]), mean_both_Dim$Wavelength, (mean_both_Dim[,2]+mean_both_Dim[,3]),
       length=0.05, angle=90, code=3, col = 'black', lwd = 2)
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_UV_dim , lwd = 2, col = "black", lty = "dashed")
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_green_dim , lwd = 2, col = "black", lty = "dotted")
lines(330:660, predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_UV_dim +  predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660)))* scaling_factor_green_dim, 
      lwd = 2, col = "black", lty = 'solid')
text(x = 650, y = 1.1, 'n = 15', cex = 1.1)
text(x = 327, y = 1.1, 'C', cex = 1.5, font = 2)

#title(main = "Fitted model on the relative sensitivity (mean±Sd) of S. noctilio insects (7 males and 8 Females) 
#      after dark adaptation",cex.main=1.5)
#legend("topright", legend = c('Fitted UV opsin model', 'Fitted green opsin model', "Sum of UV + green opsins models"),cex=1.5, lty = c(2,2,2), lwd = 2, col = c('purple', "green", 'black'))

plot(x= c(330:660), predict(model_UV_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "black", lty = "dashed", type = 'l',
     xlab = c('Wavelength (nm)'), ylab = c("Relative sensitivity (%)"), cex.lab = 0.9)
lines(330:660, predict(model_green_alpha_both_all_spectrum_2parameter_fix, newdata = data.frame(Wavelength = c(330:660))) , lwd = 2, col = "black", lty = "dotted")
text(x = 327, y = 1, 'D', cex = 1.5, font = 2)
legend('topright', c(expression(paste(lambda, 'max = 363 ± 3 nm')),expression(paste(lambda, 'max = 527 ± 2 nm'))), lwd = c(2,2), lty = c("dashed", "dotted"), col = c("black","black"))



